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/"•j ", We present a general formalism to treat slowly rotating general relativistic superfluid neutron 

stars. As a first approximation, their matter content can be described in terms of a two-fluid model, 
where one fluid is the neutron superfluid, which is believed to exist in the core and inner crust 
of mature neutron stars, and the other fluid represents a conglomerate of all other constituents 
(crust nuclei, protons, electrons, etc.). We obtain a system of equations, good to second-order in 
the rotational velocities, that determines the metric and the matter variables, irrespective of the 

X_/\ • equation of state for the two fluids. In particular, allowance is made for the so-called entrainment 

effect, whereby the momentum of one constituent (e.g. the neutrons) carries along part of the mass 
of the other constituent. As an illustration of the developed framework, we consider a simplified 
equation of state for which the two fluids are described by different polytropes. We determine 
numerically the effects of the two fluids on the rotational frame-dragging, the induced changes in 
the neutron and proton densities and the inertial mass, as well as the change in shape of the star. 

^^ ' We further discuss issues regarding conservation of the two baryon numbers, the mass-shedding 

QO (Kepler) limit and chemical equilibrium. 
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I. INTRODUCTION 

For very practical reasons, past modeling of general relativistic neutron stars has largely relied on a one-constituent 
perfect fluid approximation to describe the stress-energy of matter. This would be appropriate if neutron stars truly 
consisted solely of neutrons in a fluid state. However, this is a drastic oversimplification. In reality a neutron star is 









composed of many different constituents and is best described as a "layer-cake" . A few seconds after it is born in a 
supernova a neutron star cools below 10 10 K. Then the outer layers (with densities below, say, 1.5 x 10 14 g/cm 3 ) will 
"freeze" into a solid crust. After a few months of further cooling, the star will reach the critical temperature at which 
the bulk of the neutrons at densities above neutron drip (~ 10 11 g/cm 3 ) become superfluid. At roughly the same time 
5—1 ' ^ ne protons (expected to make up a few percent of the core material) may become superconducting. Thus, even in 
the simplest reasonable model a neutron star contains a thin fluid ocean, the roughly 1 km deep solid crust mentioned 
above, and a core containing superfluid neutrons, (potentially superconducting) protons and electrons. The lattice of 
nuclei and electrons in the inner crust is permeated by superfluid neutrons. Yet more complex models would allow 
for the presence of exotic particles in the core (eg. hyperons), high density phase-transitions to kaon and/or pion 
condensates and/or regions with deconfined quarks H . 

Despite these many complexities the perfect fluid approximation is reasonable when all the constituents flow together 
(see, for instance, Comer and Langlois H). However, it can never accomodate situations where the superfluid neutrons 
flow "independently" of the protons. Since this is, in fact, expected to be the case in a real neutron star the standard 
approach has serious deficiencies. During the last twenty years there have been continued attempts to piece together 
a coherent picture of neutron star superfluidity, using a variety of theoretical arguments and imposing observational 
constraints based on glitch and cooling rate data. Of particular interest for our present discussion are issues regarding 
the coupling between the superfluid and the normal fluid constituents, the main question being whether the superfluid 
neutrons will co- rotate with the protons or not. The answer to this seemingly innocent question is very complicated, 
but the current consensus is that even though the two fluids are strongly coupled they can sustain a state of differential 
rotation over significant timescales. 

Alpar et al [p| have shown that the charged components in the core fluid (the protons and the electrons) will be 
strongly coupled to the crust nuclei via electromagnetic interactions. This suggests that a superfluid star can be 
described in terms of two fluids: one that represents the superfluid neutrons and one that describes all charged com- 
ponents (protons, electrons, crust nuclei etc.). In the following we will refer to these two fluids as the "neutrons" and 
the "protons" , respectively. Furthermore, in a rotating star the normal fluid constituent will settle into an essentially 
non-dissipative configuration on the viscous timcscalc. This implies that the "protons" must rotate uniformly (with 



proton angular velocity £l p = constant). This is not a priori the case for the neutrons since a (pure) superfluid is 
non-dissipative. In other words, one can conceivably have a stationary configuration with differentially rotating neu- 
trons (with neutron angular velocity f2„ 7^ constant) . This would certainly be true if the two fluids were completely 
uncoupled but, in reality, this is likely not the case. In order to begin to understand what really happens we need 
to discuss a non-dissipative mechanism known as the entrainment effect and how it leads to a dissipative mechanism 
called mutual friction. 

In a mixture of the two superfluids Helium three and Helium four, it is known that the momentum of one of the 
constituents carries along (or entrains) some of the mass of the other constituent |||| . The analog in neutron stars 
is an entrainment of some of the protons, say, by the neutrons. A superfluid is by its very nature locally irrotational 
||, so when neutron stars rotate they do so via the formation of a dense array of quantized vortices. Because of 
entrainment the flow of neutrons around these vortices will induce a flow also in a fraction of the protons. This leads 
to magnetic fields being formed around the vortices. Due to the electromagnetic attraction of the electrons to the 
protons (the timescale of which is very short) there will then be a dissipative scattering of the electrons. In the context 
of oscillating superfluid stars this dissipative mechanism is known as mutual friction. 

Much of the discussion of neutron star superfluids concerns attempts to theoretically model the observed glitches. 
The favoured model has long been based on the idea that the glitches correspond to a transfer of angular momentum 
from the superfluid that permeates the crust to the lattice. This picture is attractive since 
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where Ail/O is the observed relative change in spindown rate and AI /I is the fraction of the total moment of inertia 
contained in the crustal superfluid. Furthermore, it has been argued that the time-scale on which the core superfluid 
couples to the charged component (in our case, the time on which neutron rotation rate becomes equal to that of the 
protons) is short (on the order of 1-10 s). For example, Alpar and Sauls [J7| argue that, if P is the period of rotation 
of the bulk of the matter of the neutron star in units of seconds, then on timescales of 4 x 10 2 — 10 4 P any difference 
of rotation between the neutrons and protons in the core will be damped out because of the entrainment effect and 
mutual friction. The consequence then is that only the inner crust superfluid would be able to rotate at a different 
rate from the rest of the star on a timescale beyond at most a few minutes. 

However, there are many uncertainties in this game, e.g. the strength of the vortex pinning and the exact parameters 
for entrainment, and it is not at all clear that the standard picture is correct. In fact, recent work casts some doubt 
on its validity. To see this, we follow Langlois et al g and quantify the coupling timescale in the following way: 
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where /„ and I p are the moments of inertia of the neutrons and the protons, respectively, and c r is the "drag to lift 
ratio" that depends on the relative rotation rate. This estimate is convenient because it describes the lifetime of a 
local difference in angular velocity of the two fluids in terms of a single parameter. As is immediately clear from (|l|) 
the timescale is short if c r is roughly of order unity (when the "drag" resulting from entrainment is of the same order 
as the Magnus force that acts on the vortices), while it becomes long if c r » 1 or c r « 1. Now, recent estimates 
of the strength of the pinning of the superfluid vortices to the nuclei in the inner crust suggest that c r is small in 
this region (leading to a relaxation time of ~ 50 days) || . Additionally, a detailed investigation of the scattering of 
electrons off of proton vortex clusters (Hi in the core suggests that c r is very large there. In other words, it would 
seem plausible that the superfluid components relax towards co-rotation with the normal fluid rather slowly — both 
in the crust and deep in the core. This would indicate that differential rotation between the superfluid neutrons and 
the "protons" can be sustained on timescales of days to years. 

Clearly, this is a problem where we are far away from a final conclusion at this point in time. Fortunately, this does 
not matter for the present discussion. It is clear that even if the coupling between the two fluids is strong, differential 
rotation should prevail on a time-scale orders of magnitude larger than the dynamical one (~ ms). Further evidence 
for this may follow from observed temperatures of millisecond pulsars. We recall attempts to explain the observations 
in terms of frictional heating due to differential rotation between the superfluid and the crust |llf| . This discussion 
typically suggests that roughly 1% of the total moment of inertia (eg. the superfluid in the inner crust) must rotate 
differently from the rest of the star, in such a way that 10 -2 rad/s < |fl„ — Cl p \ < 10 rad/s. This would lead to a 
significant internal heating that would be relevant after the neutrino cooling era. 

It is also worth mentioning that one should be careful before concluding that the coupling between the two fluids 
will necessarily lead to a co-rotating system. As argued by Langlois et al || , this is unlikely to happen in cases when 
an external torque is acting on the system. Then the difference J7„ — il p will stabilize at a value that represents a 



balance between the acting torque and angular momentum redistribution in the fluid. Since the external torque acts 
on the crust which then exerts a drag on the superfluid, we would typically expect fl n > fl p in an isolated neutron star 
that spins down due to magnetic dipole braking, while tt n < £l p in a neutron star accreting material from a binary 
companion. In other words, one would expect the superfluid neutron to spin faster than the charged components 
in an isolated radiopulsar that is spinning down due to magnetic dipole braking, while the opposite is the case for 
accreting neutron stars in (say) the Low-Mass X-ray binaries. Given that the latter may actually generate detectable 
gravitational waves, see |12 H], this is an interesting idea. In fact, it may well turn out that internal differential 
rotation needs to be accounted for in the construction of accurate theoretical templates for the gravitational waves 
from (say) Sco-Xl. 

During the past decades an accurate machinery for determining rapidly rotating neutron star models in general 



relativity has been developed, see 15 lq] for reviews. However, as should be clear from the above discussion, truly 
realistic neutron star calculations must account for features associated with the interior superfluid. In particular, 
realistic models for rotating neutron stars must allow the two fluids to rotate at different rates (see [Il7| for a study 
of this problem in the Newtonian context). This requires a potentially very complex extension of all existing work. 
Fortunately, given the recent development of a relativistic formalism for describing superfluids and the accurate 
numerical codes used to construct rotating single-fluid stars, many of the tools required to approach this problem 
are in hand. Hence, it is timely to initiate a program aimed at improving our understanding of rotating relativistic 
superfluid stars. 

As a first step towards the solution of this problem, the present paper concerns a formalism for modeling slowly 
rotating superfluid neutron stars. Our main aim here is to develop the mathematical framework and explore how the 
extra degrees of freedom associated with the superfluid affect slowly rotating neutron star configurations. The derived 
formalism can then serve as the starting point for relativistic studies of pulsar glitches, or spin-up of neutron stars 
due to accretion of matter. Furthermore, the present work provides the first step towards analyzing the oscillations 
of slowly rotating relativistic superfluid neutron stars. In this sense, the present work is the next logical step in the 
program initiated by Comer et al (]1| . Eventually, we aim to carry over to general relativity the pioneering calculations 
of Lindblom and Mendell Jlfl-Glll that were done for Newtonian stars. 



There are several reasons why this a ver y im portant target. First of all, several studies of the oscillations of non- 



rotating superfluid neutron stars 22 2ljj2^,|l8[ have shown that the new degrees of freedom that result because the 
neutrons flow independently of the protons lead to new sets of oscillation modes (in addition to the familiar f, p and 
g-modes) . One interesting question for follow-up studies of the present work concerns the effects of rotation on these 
superfluid modes. Another motivation for studying oscillations of slowly rotating superfluid configurations is provided 
by the gravitational- wave driven instability of the so-called r- modes |24|,|25j . A study of the r-modes in a relativistic 
superfluid star would be a natural generalisation of recent work of Lockitch, Andersson and Friedman |2q ] and would 
provide further insights into the details of such modes in realistic neutron stars. This is a crucial issue given the 
possibility that the gravitational waves from the r-modes are potentially detectable by LIGO II [ pTJJTI) . One can 
speculate that the new degrees of freedom associated with the superfluid may lead to new sets of r-modes (especially 
if the neutrons and the protons are not corotating). But even if this possibility is not realized, it is important to 
pursue an analysis of the r-modes because of the pot entially crucial dissipation due to mutual friction. This issue has 
so far only been addressed in Newtonian gravity [E8|, and there is an obvious need for a detailed relativistic analysis. 

As already mentioned above, we describe a superfluid neutron star as a two-fluid system. Our equations include 
vortices implicitly (i.e. the vorticity of the two fluids are not zero), and as argued by Comer et al |li| should 
approximate quite well the rotational properties. What are neglected are such physical effects as tension along the 
vortices, and modes that travel along individual vortices or are associated with the entire array (e.g. Tkachenko waves 
H ) . We do not assume a priori that the neutrons rotate at the same rate as the protons in deriving our equations, nor 
will chemical equilibrium be assumed. This means that we can clearly distinguish how the individual rotation rates 
of the neutron and protons (which represent independent degrees of freedom), and chemical equilibrium affect the 
rotational equilibrium. Perhaps most importantly, we wish to gain insight into how entrainment affects the structure 
(and, eventually, the dynamics) of neutron stars. This clearly requires different rotation rates for the neutrons and 
the protons since entrainment is no longer relevant when the two fluids flow together. 

Our results are presented in the following fashion: In Section II some background material will be presented 
that determines the most general, exact form for the metric and matter variables that is allowed by the underlying 
symmetries and asymptotic conditions. It will be shown that the matter field equations can be solved analytically for 
the case of rigid rotations of both the neutron and proton fluids. Section III introduces the slow-rotation approximation 
and applies it to the superfluid and Einstein field equations. In Section IV we discuss the numerical solutions to the 
field equations, and discuss results for a simple model equation of state. Finally, we discuss issues regarding "physical" 
sequences of stars, such as baryon number conservation, the Kepler limit and chemical equilibrium in Section V. Our 
concluding remarks are offered in Section VI. Some further details that are useful to the derivations discussed in the 
main text are presented in Appendix I. Except where explicitly noted we use units such that G = c = 1. 



II. AXISYMMETRIC, STATIONARY, AND ASYMPTOTICALLY FLAT SPACETIMES 

The study of axisymmctric, stationary, and asymptotically Sat spacetimes has by now a long history. This means 
that there is a considerable literature on the subject, and we refer the reader to [ |15|Jl6| for detailed reviews. Our modest 
constribution in this paper is to construct rotating superfluid stellar models within the slow-rotation approximation. 
Key references are Comer et al Cq] and Langlois et al M for the superfluid formalism, Bonazzola et al p9] for how to 



set up axisymmetric, stationary, and asymptotically flat spacetimes, and the work of Hartle and collaborators []30 31 
for the slow-rotation approximation in relativity. 

A. The Conditions of Axisymmetry, Stationarity, and Asymptotic Flatness 

Axisymmetric, stationary, and asymptotically flat spacetimes are defined in the following way fl29[ : (i) There exists 
a Killing vector, to be denoted here as t p , that is timelike at spatial infinity; (ii) there exists a Killing vector, to be 
denoted here as ^ , that vanishes on a timelike two-surface (called the axis of symmetry), is spacelike everywhere 
else, and whose orbits are closed curves; and (iii) asymptotic flatness means the scalar products t v t v ', <$> v <$> v and t u (jf 
go to, respectively, — 1, +00, and at spatial infinity. We will likewise impose the so-called "circularity condition" on 
the stress-energy tensor T p , which means that 

TX = oct" + (3^ , Tift" = At" + g^ . (2) 

Carter |p2| has shown that assumptions (i), (ii), and (iii) above imply that the Killing vectors commute. It thus 
follows that a coordinate system, x° = t, x , x 2 , and a; 3 = <fi, exists for which 

t" = (1,0,0,0) , 0"= (0,0, 0,1). (3) 

Carter p3] has also shown that the "circularity condition" implies that the metric tensor can be reduced to 



g^udx^dx" = -N 2 dt 2 + g u (d</> - N^dtf + g n (da; 1 ) 2 + 2g 12 dx 1 dx 
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where each of the components depends on only the coordinates a; 1 and x 2 . By choosing so-called orthogonal coordi- 
nates, one may further simplify the metric so that g±2 — without any loss of generality. 

There are a few different orthogonal coordinate systems that have been used to study the types of spacetimes under 
discussion here. Bonazolla et al |2£J discuss several examples that have been used for exact calculations, and for 
studies using the slow-rotation approximation. For the present discussion, we will use orthogonal coordinates x 1 = f 
and x 2 = 8 common to slow-rotation studies so that the metric is 

g^dx p dx» = - (n 2 - sm 2 6K [iV ] 2 ) di 2 + Vdf 2 - 2sm 2 6KN 4 'dt d<p 

+K(d0 2 +sin 2 9d<f) 2 ) . (5) 



B. Axisymmetric, Stationary, and Asymptotically Flat General Relativistic Superfluid Neutron Stars 

We want to see what the various conditions discussed above imply for a superfluid neutron star. Thus we first recall 
the formalism that has been used to model general relativistic superfluid neutron stars |l|,0] . The central quantity is 
the so-called "master" function A, which is a function of the three scalars n 2 = —n p n p , jr = —p p p p , and a; 2 = —p p n p 
that are formed from n M , the conserved neutron number density current, and p p , the conserved proton number density 
current. The master function is such that —A corresponds to the total thermodynamic energy density. 

A general variation (that keeps the spacetime metric fixed) of A(n 2 ,p 2 , a; 2 ) with respect to the independent vectors 
n p and p M takes the form 

5k = n P 8n" + Xp Sp p , (6) 

where 



Hfj, = Bn^ + Ap^ , Xm = C £V + -4»fy > ( 7 ) 

and 

The covectors /i^ and x M are dynamically, and thermodynamically, conjugate to n^ and p 11 , and their magnitudes 
are, respectively, the chemical potentials of the neutrons and the protons. The two covectors also make manifest the 
entrainment effect since it is seen explicitly how the momentum of one constituent carries along some mass current of 
the other constituent (for example, /i^ is a linear combination of n^ and p^). We also see that there is no entrainment 
if the master function A is independent of x 2 (because then .4 = 0). 
The stress-energy tensor is given by 

rf = f«+^ + »V, (9) 

where the generalized pressure \& is given by 

* = A - n>„ - p p Xp ■ (10) 

The equations of motion consist of two conservation equations, 

V M n" = , V^ = 0, (11) 

and two Euler type equations, which can be conveniently written in the compact form 

^V [M /v ]= , p"V[^]=0. (12) 

When all four equations are satisfied then it is automatically true that V M T^ = 0. 

It is convenient for what follows to rewrite each of n M and p^ as a product of a magnitude with a unit timelike 
vector in such a way that 

n^ = nu^ , pf" = pvP , (13) 

where u^u^ = — 1 and v^v^ = — 1. One can easily show that the circularity condition, when applied to the stress- 
energy tensor in (O), leads to the following forms for the unit timelike vectors u^ and v^ (using the form of the metric 
given in (||)): 

f* + <W „, M _ t» + n P r (14) 



Jn 2 - sm 2 9K (AT* - n n f ^N 2 - sm 2 6K (JV* - Vt p f 

where Q n and £l p are the angular velocities of the neutrons and protons, respectively. 

With the appropriate decomposition of the unit vectors established, the remaining matter variables can be simil arly 
decomposed, and the matter field equations can be analyzed. It turns out that the two conservation equations (|ll| ) 
are automatically satisfied. The other matter field equations are the two Euler relations dT2) , which can be shown to 
reduce to 
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sin0\/F (N* - fi n ) sin6»VX (TV* - Q p ) 
"n= ^ , lo p ^ . (17) 

We will focus most of our attention on the case of rigid rotation. Then Q, n and Vt p are both constants and the Euler 
equations lead to the following first integrals of the motion for the neutron and protons, respectively: 

Ap [1 - LUn^p] 



N Bn^/T~ t 



Xc = N [c P f^ p + An ^-_^A , (18) 

where fi c and \c are both constants. When the rotation is set to zero for both fluids, we retain the first integrals 
obtained by Comer et al |lg| |. 

Given the above results the superfluid field equations have been completely solved. However, we still need to 
determine the corresponding metric functions. To do this we must consider the Einstein field equations. Unfortunately, 
these are sufficiently complicated that the same level of progress in solving them can not be achieved. That is, we 
cannot write down the solution in closed form. Such a solution was obviously not expected since the Einstein equation 
requires a numerical solution already in the one-fluid case. 

The equations we have discussed so far do not include any approximations (apart from the particular two-fluid 
model for the superfluid). This means that they will be relevant irrespective of the stars rotation rate. Thus, they 
could in principle be used as a basis for a numerical solution for rapidly rotating stars, following in the footsteps of p9[ . 
However, we want to explore the new degrees of freedom that come into play when we consider a two-fluid system. 
To do this, it seems natural to first consider the problem in the slow-rotation approximation where the equations are 
somewhat more transparent and we can make further "analytical" progress. We will return to the problem of rapidly 
spinning stars in the future. 

In the following Section we will introduce the slow-rotation approximation and derive the relevant field equations. 
These will subsequently be integrated numerically, with sample results being discussed in Sections IV and V. 

III. THE SLOW-ROTATION APPROXIMATION 

In order for the slow-rotation approximation to be valid, the angular velocities must be small enough that the 
fractional changes in pressure, energy density, and gravitational field due to the rotation are all relatively small. 
When applied to our system the approximation can be translated into the inequalities, cf. |30f| 

il 2 n or 0* or n n n p « (-^) 2 ^ , (19) 

where the speed of light c and Newton's constant G have been restored, and R and M are the radius and mass, 
respectively, of the non-rotating configuration. Since GM/c?R < 1, the inequalities also imply 

Q n R « c and Q p R « c . (20) 

These conditions indicate that the slow-rotation approximation ought to be useful for most astrophysical neutron 
stars. 

For example, we can compare the above conditions, e.g. that 

[n n> n p> V*U^«y^r« 11500 (^) (jnp) s_1 ' (21) 

to the empirical estimate for the Kepler frequency (i.e. the rotation rate at which mass-shedding sets in at the equator) 
that has been deduced from calculations using realistic supranuclear equations of state Jlq] : 
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We should also recall that the fastest observed pulsar rotates with a period of 1.56 ms. This corresponds to £1 rs 
4000 s _1 , i.e. roughly half the Kepler rate. In view of these examples, it is not too surprising that calculations for 
one-fluid models have shown the slow-rotation approximation to be accurate to within (say) ten percent for stars 
spinning at the mass-shedding limit. 

Before combining the slow-rotation approximation with our superfluid formalism, it is worthwhile emphasizing two 
key points that follow since we are dealing with two distinct fluids: (i) If only one of the rotational directions is 
reversed, then the physical configuration of the star should change; but (ii) if both rotational directions are reversed, 
then the physical configuration should not change. As it turns out, the only quantities that contain terms linear 
in the angular velocities are the metric coefficient N^ , that represents the dragging of inertial frames, and the fluid 
four- velocities u^ and v^. Thus all other effects due to rotation enter only at the second-order in the angular velocities. 
The calculations performed here, then, will keep terms through second-order. 

A. Slow Rotation Expansion of the Metric and Matter Variables 

Not very surprisingly, the equations that determine the metric variables in the slow-rotation approximation for our 
two-fluid system are similar to the one-fluid ones derived by Hartle more than 30 years ago p0[ . Of course, we want 
to allow the two fluids to rotate at different rates so there are some conceptual differences associated with the fluid 
variables. Still, it is natural to follow Hartle |3Q] and expand the five independent metric coefficients as 

N = e"^ /2 (l + h(r,6)) , 

V = e A(f) (1 + 2v(r , 0)) , 

K = f 2 (l + 2fc(f,<9)) , 

N+=w(r,9), (23) 

where it is to be understood that the terms h, v and k are each second-order in the angular velocites whereas the 
frame-dragging u) is a first-order quantity. Initially we will assume that the two fluids rotate rigidly, which means 
that there is no need to write explicitly similar expansions for the fluid velocities. Thus, for the fluid we need only 
define the slow-rotation expansion for the neutron and proton number densities n and p, respectively: 

n = n (r)(l + r)(r,6)) , p = Po (f) (1 + fc(f,0)) , (24) 

where the terms r\ and $ are understood to be second-order in the angular velocities and we have introduced the 
convention (to be used throughout the paper) that terms with an "o" subscript are either contributions from the 
non-rotating background or quantities that are evaluated on the non-rotating background (e.g. x\ — n p etc.). The 
expansions of the remaining fluid variables (such as the velocities, the stress-energy tensor components, etc.) can all 
be obtained in terms of the metric and particle number density relationships written above (some of the more useful 
results are presented in Appendix I). 

A non-rotating background configuration is specified once n , p Q , v and A are known. The two background metric 
components are obtained as solutions to 

A' = ^—J- - 8nre x A , v' = -^ — + 87rre A * (25) 



while the neutron and proton number densities are obtained from (see |18 for a complete discussion 
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where 
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Regularity of the geometry at the center of the star requires that A(0), A'(0), ^'(0), n' o (0), and p' o (0) all vanish. The 
surface of the star is the value R of the radial variable for which ^ Q (R) = 0. Finally, the total mass M of the 
configuration is given by 

r R 
M = -4tt / f 2 A (f)df . (28) 
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The equations that determine the rotational features are obtained by taking the expansions given above and putting 
them into the full superfluid and Einstein field equations, but keeping only terms up to second-order in the rotational 
velocities. Even with the slow-rotation approximation, this new set of equations represents a two-dimensional problem. 
Fortunately, the dimension represented by the 9 coordinate can be successfully taken care of by an expansion in 
Legendre polynomials Pi(cos9), and then it can be shown that we need only consider the I = 0, 1,2 components (the 
argument is completely analogous to that of Hartle [[30|). The net result is that the metric corrections h, v, and k can 
be written as 

h = h (r) + h 2 (r)P 2 (cos6) , 

v = Vo(f) + v 2 (r)P 2 (cos9) , 

k = k (r) + k 2 {?)P 2 (cos9) , (29) 

where P2(cos#) = (3cos 2 6>— l)/2. Futhermore, Hartle |3(| has shown that a coordinate transformation can be imposed 
whose sole effect is that fco(f) vanishes. As for the remaining metric piece u>, it can be shown that it is independent 
of 9, exactly as in the one fluid case, so that we can write ui = ui(r). It will be useful for what follows to define the 
quantities 

L n = u> — Q n , L p = lj — ftp . (30) 

Up to an overall minus sign, these represent the two rotation frequencies as measured by a local zero-angular momen- 
tum observer. 

The remaining two matter fields 77 and $ are written as 

V=Vo(r)+m(r)P2(cose) , $ = * (r) + * 2 (r)P 2 (cos6») . (31) 

It is also convenient at this point to write each constant that appears in the first integrals of the Euler equations p2) 
as a sum of two terms, one (e.g. /Xqo) which is for the non-rotating background and a correction (essentially 7„) which 
is second-order in the rotational velocity, i.e. 

Mc = A*oc (1 + In) , Xc = Xoo (1 + lp) ■ (32) 

To complete the system of equations that determines a slowly rotating superfluid model, we need the Einstein 
field equations. The field equations that will be used here are G% = 8vrT °, G\ = 871^, i? 3 = 8tt (T 03 - [l/2]T ff0 3), 
^2 — -^3 = ^ (T 2 — T3), and i?i2 = 0. The only matter field equations that need to be considered are the first 
integrals in (Eq). It should be noted that not all of these field equations are independent; a proper subset will be 
extracted in Sec. IV where we discuss the numerical solutions. 

Finally, we need to address a subtle point (discussed in detail by Hartle 130]): Near the surface of the star the 
conditions will be such that, for instance, the ratio of perturbed to non-perturbed neutron number density will be 
ill-defined because the non-perturbed background density will vanish. That this will be the case is easy to understand 
since the centrifugal force changes the shape of the star. To overcome this complication, we introduce a new radial 
coordinate r given by 

r = r + £(r, 9) (33) 



such that 



A(f(r,0),0)=A o (r) 



(34) 



That is, r measures the average distance from the centre to the rotational surfaces of constant energy density [J34J 
The situation is illustrated schematically in Figure nL We then expand the displacement function £ as 



£ = £o(r)+£ 2 (r)P 2 (cos0). 
The shape of the surface of the star can now be obtained as 

R(9) =R+£(R,6) . 
To linear order in £ we now find that 

A(r,0)=A o (r)-A o (r)£(r,0), 



(35) 
(36) 



(37) 



where the prime denotes a derivative with respect to r. This and Eq. (105) from Appendix I imply for I — that 

,•2 



while for I — 2 we have 



V 9 

fJ- n r) Q + XoPo^o + tt^ A n p (fi„ - Q p ) = A^ 



T 2 

^ n Tl2 + XoPo®2 - T^AnoPo (H n - ^p) = K& 



(38) 



(39) 




FIG. 1. A schematic illustration of the rotationally distorted constant density surfaces and the meaning of the coordinate 
transformation f — > r + £ that is used in the slow-rotation case. 



The coordinate transformation affects the "03", "22 - 33", and "12" Einstein equations and the first integrals of 
the Euler equations only to the extent that f gets replaced with r. More substantial is what happens to the "00" and 
"11" Einstein equations. For these one must transform the corrections to the Einstein and stress-energy tensors in 
the old coordinate system, denoted by <5G£J and 5T^, respectively, to the corrections in the new coordinate system, 
denoted by AG£J and AT^ . The transformations are readily found to be 



AGg = «G$(r,0) + ^W(r,0) 



Or 



for the Einstein tensor corrections and 



AT ° = 5T°(r,0) 



dr 



H(r,0) 



AG{ =5G{(r,e) 



ATI =STl{r,9) 



dG\\ 

1 lo 



Or 



£M) 



dTl\ 

1 lo 



Or 



Z(r,6) 



for the stress-energy tensor corrections. 

Now, the first integrals of the Euler equations have as their I = contributions 



ROl „ a0\ n 



Mo 



r 2 p ( . dA 



+ n p 



dA 
dx 2 



yiin ^"p) 



(40) 



(41) 



3e< 



;L% + h Q 



C °\oP°^ A ol n ° r 2 n / 

7 P = -^^*o + '° Vo + t^j— [Ao+po 
Xo Xo 3e» Xo V 



dA 
dp 



n p 



dA 
dx 2 



\^"Tl t^p) 



where fi = B n + A p and Xo = C p + A p , and for I = 2 

ff °| o no ; 4L^ a r 2 Po ( dA 

= — m H 2 — $ 2 - TT — [Ao + n — 

Mo Mo 3e" Mo V on 



™oPo 



a4 



^i£ n lip J 



(42) 



3^™ + ^ 

Xo Xo 36^ Xo V "P 



™oPo 



&4 

dx 2 



\Mn Mp) 



^ + /l2 



The "03" Einstein equation yields 

^ (r 4 e -( A +")/ 2 Z; i )' - 167re( A -^/ 2 (*„ - Ao) L n = 167r e ( A -^ 2 XoPo (fi„ - Op) 



(43) 



(44) 



which is the equation that determines the frame-dragging. It is formally equal to that obtained by Hartle |3(J except 
for the non-zero source term on the right-hand-side. In the particular case when the two fluids are co-rotating 
(fi„ = ftp) we retain the standard result. 

The "22 - 33" and "12" equations only have I — 2 contributions and these are, respectively, 



v 2 + h 2 = 



Q e U+X 

87rr 4 
+ 3e" 



I' \ 2 , 87rf4 at, A n ?a 



,, ( , (*o-A )^ 

XoPo (fin - fip) ( in + £p ) - A>™oPo (fin ~ fip) 



and 



- (% + M - (*a + h 2 )' - V — (h 2 -v 2 )=Q 
r 2 



The I = part of the "00" equation is 



I6nr 2 



(*o - A ) l£ + XoPo (fin - fip) [Ln + L p \ - A Tl p (fi„ - fip)^ 



8-AoCo-^(^o)' + ^ T (i;. 



with the I — 2 piece being 
lQirr 2 



= 



3e" 



(*o - A ) L 2 n + XoPo (fin ~ fip) [L n + Lp) - A n p (fi n - fi p )' 



2/7* \ ' r / ~ \ * 

8 ^ A ^ 2 - ^ b W2 ) - ^TX fa) + 

Finally, the I — part of the "11" equation is 



2 2 

7\ 



fc 2 



A' 






2e> 



-,v 2 . 



(45) 



(46) 



(47) 



(48) 
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2 2 / , 1 

\ x 



-+-j-o + ^ TI 



^n ) = 8?r [Mo^o7n + XoPolp ~ (^o ~ A Q ) fc() + 



— [jj, n L n + XoPoLlj - — n p ./4o (H„ - Op)' 



(49) 



while the Z = 2 result is 



1-fc' - - 

re r 



2 
re > 



— 5-/l' 2 - -jfca r(l/ + -)«2 + - T (l> / + -)*2- ""ofe ~ 



4 



6e i/+A 



(*)' 



2 

(* - A ) h 2 + — {jJ,on L 2 n + XoPoLlJ 



3e^ 



(50) 



Before we proceed, it is worth emphasising the fact that the I = 0, 1, 2 components of the fluid and metric field 
variables decouple. This is in complete analogy with the results of Hartle for the one-fluid case |30(] . 

B. The Exterior Vacuum Solutions 

In the exterior of the star our problem is identical to the one-fluid one. Hence, we can use the analytic solutions 
found by Hartle EG]. Thus, the vacuum solutions are written 



2J 



(51) 



ho{r) = — + 



r-2M r 3 (r-2M) ' 



(52) 



vo(r) = -ho(r) = 



J 2 



SM 



r-2M r 3 {r-2M) 



(53) 



Mr) = -A 



3,rx>/ 2M\ r_2M 

2 \mJ \ r I \ r 



(r - M) (3 - 6M/r - 2{M/r) 2 ) 
M (1 - 2M/r) 



Mr 3 \ r J 



(54) 



k 2 (r)=A 



21m 1 U " 



2AP 



In 1 - 



2M\ 3{r- M)~8{M/rf(r- M/2) 



M (1 - 2M/r) 



J 2 



Mr 3 



2M 
r 



(55) 



v 2 (r) = A 



3 / r \ 2 / 2M\ , / 2M 
2(m) ^- W *-- 



(r - M) (3 - 6M/r - 2(M/r) 2 ) 
M(l - 2M/r) 



J 2 / 5M 

Afr 3 l r 



(56) 



where J, SM and A are constants. J represents the total angular momentum of the star with respect to inertial 
observers at spatial infinity while SM represents the change in total mass-energy due to the rotation. The final 
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constant A is related to the quadrupole moment of the star through (cf . Hartle and Thorne |31j , although with the 
sign reversed to agree with Laarakkers and Poisson |35f| ) 

and it is determined when we solve the I = 2 equations in section IV. This then facilitates a comparison with the 
well-known result for the Kerr black hole spacetime. Namely that 

QKerr = ~ ■ (58) 

Assuming that the metric is continuous at the surface of the star, it can be shown that 

J = -Yf r 4 e {X - )/2 (*o - A ) L n dr - ^ /"* rV^>/ 2 XoPo (O n - Q p ) dr = 



~ [ R drr\' X -'■' 



ii 





fi n L n + XoPoL p . (59) 



It would be tempting to interpret the two terms in the integrand as representing the angular momentum in the 
neutrons and the protons, respectively. However, as we will show in Section VB this is only correct in the absence of 
entrainment. We also have 

I 2 
SM = (R- 2M)v (R) + -3 . (60) 

Using these two equations, the dependence of J and SM on the central neutron and proton number densities can be 
established. 

IV. NUMERICAL RESULTS 

There are four essential steps to using the slow-rotation formalism to produce a model of a general relativistic 
superfluid neutron star: (i) build a nonrotating background configuration, (ii) determine the frame-dragging by 
solving ( J44| ) using as input the background configuration, (iii) solve the I = equations using as input the background 
configuration and the frame-dragging results, and (iv) determine the I = 2 contributions by solving the relevant 
equations in a similar way. In this section we describe our implementation of these steps and present some typical 
numerical results. 

A. The Background Configurations 

The nonrotating background configuration is easily constructed by solving equations ( p5| ) and (Eq). The required 
parameters are the central number densities n o (0) and j> o (0). For any given equation of state these parameters are 
linked by the requirement that the background fluid be in chemical equilibrium. As a suitably simple model equation 
of state, we consider the case when the two fluids are described by independent polytropes. The master function then 
takes the form |18| 

A(n 2 ,p 2 ) = -m n n - cr„rA - m n p - a p p Pp (61) 

For simplicity we assume that the neutron and the proton masses are identical (— m n ). It should also be noted that 
we do not, from this point on, include the entrainment effect, even though all equations we have derived so far, in 
principle, accomodates it. The main reason for not building specific numerical solutions with entrainment here is 
that the mere inclusion of two fluids is significant and we feel it is crucial to first establish that the slow-rotation 
framework that we have developed produces reliable results. Moreover, the implementation of entrainment requires a 
considerable amount of further work as far as the model equation of state is concerned. In practice, the fact that we 
are not including entrainment in the models discussed in this and the following section means that all our numerical 
calculations assume that A = Aq = etc. A future paper will be devoted to a model for entrainment and its effect 
on rotating neutron stars. 

12 



Our model equation of state is the same as that used by Comer et al jT|] in the study of quasinormal modes of 
non-rotating superfluid stars. However, we have chosen our parameters differently in order to make our star somewhat 
more realistic. Our aim was to create a model star with mass roughly 1.4M , radius 10 km and a central proton 
fraction of roughly 10%, i.e. what could be considered typical values for a realistic neutron star. Given such a model 
we can compare the numerical results for (say) the increase in inertial mass due to the rotation with results obtained 
for more realistic supranuclear equations of state and stars with realistic parameters. 

The particular values we use for the two polytropes are 

O n = 0.2 , 

n = 2.3 , 

o p = 2 , 
[3 P = 1.95 . 

We assume that the non-rotating model is in chemical equilibrium. This means that we impose the condition 
Mo = Xo throughout the star. The resultant sequence of stellar models is illustrated in figure o. The particular 
model which will be used to illustrate the typical effects of rotation corresponds to a central neutron number density 
n o (0) = 0.93 fm~ 3 , which leads to a total central energy density of — A o (0) = 1.215 fm -3 . The resultant non-rotating 
star has M = 1.409Af Q and R — 10.076 km. As is clear from figure || this star is stable against radial oscillations. 
The neutron and proton density profiles for our particular model star are shown in figure 0. 
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FIG. 2. The neutron and proton number densities as functions of r for our model star. 

This figure illustrates two important things. First of all it is easy to see that we have chosen our parameters in 
such a way that the proton fraction at the centre of the star is roughly 10%. Secondly, it should be noted that the 
proton and neutron densities vanish at the same value of r. In other words, the surface of the star corresponds to 
n (R) = Po(R) = 0. This is, in fact, a consequence of our choice to enforce chemical equilibrium throughout the star. 
It obviously means that we are not at this point considering the effects of regions of the star that extend beyond the 
superfluid: the low-density crust region and the ocean expected to cover the neutron star's surface. To work with 
a simplified stellar model is natural at this point, but it is worth pointing out that a formalism for relativistic solid 
crusts was developed by Carter and Quintana [ p7[ a long time ago. A potentially very interesting extension of the 
present work concerns the matching of our equations, representing the superfluid core, to a slow-rotation version of 
the equations of Priou |3£| , describing the dynamics of solid crusts. 

B. The Frame-Dragging 

Once the background configuration has been determined, the frame-dragging lo can be calculated from (MJ). This 
calculation is standard, but it is notable that the two-fluid results differ from those of Hartlc in that a single integration 
does not suffice to determine the frame-dragging at all rates of rotation. This is simply because the integration of 
( [44| ) requires the central value of L n as well as both the neutron and the proton rotation rates. In practice, we solve a 
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rescaled version of ( EJJ ) such that we determine L n = L n /£l p given L n (0) and the relative rotation rate £l n /Q p . This 
then determines the frame-dragging for all £7 p with a fixed relative rotation. It is natural to scale the variables with 
the rotation rate of the protons rather than that of the neutrons since, as we discussed in the Introduction, one would 
expect the charged components in a neutron star core to be electromagnetically coupled to the nuclei in the crust. 
This locks the two components (which make up our "proton" fluid) together on a very short timescale. Meanwhile the 
superfluid neutrons may rotate differently. Furthermore, as will become clear in Section V, there may be situations 
where one would need to allow the neutrons to rotate with f2 n ^ constant. 

The solution to the frame-dragging equation (|4J) is to be such that the interior matches smoothly onto the known 
vacuum solution (plft. This means that we must have 



L n {R) = — fl n 



■2J_ 



(62) 



We can easily see that L n and its derivative are smooth provided that we have 



L„{R) — — fi„ — — — : — 
3 dr 



(63) 



r=R 



Once we find a value for L n (0) such that (B3h is satisfied, we have an acceptable solution to the problem, and we can 
determine the angular momentum of the configuration from (|62j) . 

Having solved the frame-dragging equation for our model star we find (not very surprisingly) that the result differs 
very little from the standard one-fluid result as long as the relative rotation of the neutrons and the protons is not 
too extreme. One typically finds that the frame dragging is a smooth function that decreases monotonically from the 
centre of the star towards the surface. A sample of results are shown in figure ||. 
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FIG. 3. Typical results for the frame-dragging ui for superfluid stars in which the neutrons and the protons rotate at slightly 
different rates. 

Even though the outcome may be of limited physical significance it is interesting to consider more extreme situations, 
such as ones where the neutrons are rapidly counter-rotating relative to the protons. An example of such a case is 
shown in figure 0. Here the neutrons rotate backwards with respect to the protons in such a way that Q n = — 0.08Sl p . 
This example shows that the two-fluid system allows models where the frame-dragging no longer changes monotonically 
from the centre to the surface of the star. In the particular case shown in figure the forward rotation of the protons 
(which make up roughly 10% of the total number density at the centre, cf. figure 0) counteract the frame dragging 
induced by the neutrons in the central parts of the star. Meanwhile, the backwards rotation of the neutrons dominate 
the frame-dragging in the surface region. 
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FIG. 4. The frame dragging in an extreme (rather unphysical) situation where the neutrons and the protons counterrotate 
in such a way that the net frame dragging is "backwards" at the surface of the star but "forwards" in the central parts. 



C. The / = Results 



Let us now turn to the next phase of the problem, namely the I = equations at the 0(f2 2 ) level. Thus we consider 
Eqns. (|38|), (42), ( (47|) and (49). This corresponds to five equations for £ > ?7o> &o, h and vq- We also need to determine 
the two integration constants 7„ and j p . In practice we proceed as follows: First we introduce a new variable 



mo = re vq 



(64) 



Then we implement the boundary conditions at the centre of the star by requiring that mo(0) = £o(0) = and 
ho(0) — constant. The first two correspond to spacetime being locally flat at the centre and the requirement that our 
new coordinate system has the same origin as the spherical coordinates used to describe the non-rotating model. It is 
also worth noticing that the condition £o (0) = leads to the central energy density of our rotating star necessarily being 
equal to — A o (0), i.e. the central density of the non-rotating star. This is true also for the single fluid slow-rotation 
approximation, cf. [p0[ . 

However, in the superfluid case we have an extra degree of freedom. It is not sufficient to provide the total energy 
density at the centre, we need to also provide either t/q (0) or $0(0). These represent the change in the neutron and 
proton central density, respectively. The fact that we have £o(0) = relates these two parameters and we must have 



(Von o r) )\ r=0 + (XoPo$o)| r=0 = 



(65) 



In other words, we only need to provide rjo(0) (say). In the following we will typically assume that ?7o(0) = <&o(0) = 0, 
i.e. that the central proton to neutron ratio remains unchanged in the rotating star. However, this choice is made 
for convenience and there are other options that may be physically more relevant. We discuss such possibilities in 
section V, but all results presented in this section pertain to the case when ?y(0) = <&o(0) = 0. 

Given the data at the centre of the star we can easily determine j n and 7 P from ( f42|) . This leads to 



h (0) 



> .4°| o -Xoffo 



°Po*0 



7 P = h Q (0) 



' I, r°\ — v a 



M0X0 



^0*0 



(66) 



r=0 



Now we are set to solve the various I = equations. We do this by integrating (|47|) and (|49|) using a standard fourth 
order Runge-Kutta scheme. After each integration step we determine the relevant values for £o> tjq and <&o by solving 
the algebraic equations (j^) and (42). At the surface of the star the solutions for ho and mo must be matched to the 
vacuum solutions given in section IIIC. Once we find a value for ho(0) such that both ho and its derivative match the 
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vacuum result we have our desired solution to the I = problem. The matching at the surface also determines SM, 
the rotationally induced change in inertial mass of the star. For the particular stellar model shown in figure |2| and 
the case when the two fluids are corotating we find that 



SM = 0.091 



VI kl 



kHz 



Mr. 



where v v — tt p /2ir is the rotation rate of the protons. As we will discuss in Section VB the mass-shedding limit for 
this star corresponds to v v — 1630 Hz, which means that rotation may increase the inertial mass of the star by up 
to roughly 17 %. The rotationally induced increase in mass for a sequence of models is illustrated in figure |5[ We 
compare the mass of each non-rotating star to the maximum mass (which is reached when the star spins at the Kepler 
limit (see section V for a discussion). In the figure we indicate (by an arrow) the rotating analogues of the model 
star from figure that we construct within the slow-rotation approximation. We also show more "realistic" sequences 
based on conserving the total baryon number (see section V). In addition, the figure shows that, just as for single-fluid 
stars, there will be a family of supramassive stars that have no stable non-rotating counterpart J36|. 
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FIG. 5. We compare the mass of the non-rotating star (solid line) to its maximum value in the rotating case (which is 
reached when the star spins at the Kepler limit, dashed line) for a sequence of stars parameterized by the central energy density 
— A o (0). The vertical arrow indicates rotating configurations with the same central energy density as the canonical model shown 
in figure B. The dotted line starting at the base of the arrow shows a more "realistic" sequence of rotating stars constructed by 
conserving the total baryon number. We also show a similar sequence corresponding to the maximum mass non-rotating star 
(upper dotted line). Above this latter line, the various stars have no stable non-rotating counterpart, they are the so-called 
supramassive stars. 

Another illustration of the effects that rotation has on the inertial mass of the star is shown in figure q. Here we 
show how the total mass for a star rotating at the Kepler limit depends on the relative rotation between the neutrons 
and the protons. The particular results are for our canonical stellar model from figure @, and the behaviour is easily 
understood if we compare the results for the mass to the Kepler frequency results in figure \L2. 
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FIG. 6. An illustration of how the total gravitational mass for a star rotating at the Kepler limit depends on the relative 
rotation between the neutrons and the protons. These particular results are for our canonical stellar model. 

The obtained solution for ho is illustrated in Figure R. The variables £o , rjo and $o that are solved for algebraically 
are shown in figures @-[L0[ In these figures it is worth noticing that the rotational correction to the neutron number 
density 770 diverges at the surface. One can easily show that such a divergence will be present whenever the index of 
the corresponding polytrope (here /?„) is larger than 2. It should be pointed out, however, that this does not affect 
the determination of any physical parameters such as the total mass or the corresponding baryon number. Finally, 
we note that the function rao increases monotonically from the centre to the surface. Hence, we do not illustrate it 
here. 
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FIG. 7. The metric functions /io, /12 and k for our model star and the case when the superfluid neutrons corotate with the 
protons (fi n = fl p ). This particular data corresponds to the protons rotating at v v = 1 kHz. 



D. The / = 2 Results 



the I = 2 equations @, (g|), 
3f) contains no new information 



The final step in constructing a slowly rotating superfluid star consists of solvin 
(|45|), [46|) and (|5C|). These equations determine ^2, ?72, $2, ^2, V2 and h-i. Equation ( 
and we will not use it here. 

The strategy for solving these equations is similar to that we adopted for the I — problem. In the / = 2 case it 
is, however, useful to rewrite the equations somewhat. We note that once a solution to the differential equations ( |46| ) 
and (|50|) is obtained, all other variables follow algebraically. As Hartle points out |3Q], it is convenient to introduce a 
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new variable k = ki + hi . This definition and some straightforward manipulations lead to two coupled equations for 
k and ho. These can be written 



k' + u'h 2 = — (2 + rv') 




47r r 
~3 



(*o- Ao)Ll+xoPo{L n + L p ){n n -n p )- A n oPo (n n - n p ) 2 } , (67) 




+ 8^e A (A -* ) 

2 



4e A r 



- (r 2 e-V) 2 -2) 




+ 167rr 2 (zy') 2 (*o-A ) J L ? 2 l 



+ 32ire x (n n Ll + X oPoL 2 p ) + l6n X oPor 2 ^') 2 (^n - %){L n + L p ) - 16Trn oPo A (2e A + r 2 {v') 2 ) (Cl n - O p ) 2 } . (68) 



By analyzing these equations for small values of r we find that in order to have a regular solution we must require 
that \i2 — c\r 2 and k — C2r A near the centre of the star. A series expansion of ( p8| ) then reveals that the two constants 
c\ and ci should be related according to 



C 2 + 

2tt 
3^ 



2vr ($0- -Aoj ci = 
2(* - Aq)! 2 - n p A (n n - n p ) 2 ~ 2xoPo(^n - £l v )(L n + L p ) - (/io«o^ 2 + XoPoLl) 



(69) 



Furthermore, since we know that £2 — » as r — > we can deduce from (|39j) and |43|) that we must have 772 = $2 = 
at the centre of the star. 

Schematically, we solve the I = 2 equations as follows: First we find a solution to (g7|) and (|68| ) that matches 
smoothly onto the vacuum solution at the surface of the star (cf. Section IIIC) . We do this by finding a particular 
solution to the full problem and adding to it a multiple of the solution to the corresponding homogeneous problem, 
where the right hand sides of (|67| ) and (68) are set equal to zero and the two constants are related by 



c 2 + 2tt ( * - r A ) d 



(70) 



The desired solution for [k, /12] is then the linear combination that gives the vacuum relation between these two 
variables at the surface. Having determined the appropriate linear combination of the interior functions (as well as 
the constant A in the equations in Section IIIC) we can readily reconstruct the full solution, and in addition solve 
( p9| ) and (J43|) for £2, 772 and $2 after each integration step. This then completes the construction of a slowly rotating 
superfluid stellar model. Typical results for our model star are shown in figures |7|-fU)| 

Once we have solved all the / = 2 equations we can work out the rotationally induced change in the shape of the 
star, i.e. the centrifugal flattening. This effect is conveniently expressed in terms of the ratio between the polar and 
equatorial radii, R p and R e respectively. Within the slow-rotation approximation this leads to 



Rp 



36(A) 
2R 



and for our model star (with 7/0(0) = and Q n = Q p ) we get 

3e „ 1 _ 113 r~^-V 

Thus, we find R p /R e rs 0.7 for a star spinning at the mass-shedding limit. 



(71) 



(72) 
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FIG. 8. The functions £o and £2 for our model star and the case when the superfluid neutrons corotate with the protons 
(Q n = f2 p ).This particular data corresponds to the protons rotating at v v = 1 kHz. 







0.5- 


_____^J 


0.0 




0.5- 




""-.^ 


1 0- 




1 




^0 


- 


n o ri 2 


1 


1.5- 


, 


1 1 1 1 1 1 1 



0.0 



0.2 



0.4 0.6 

r/R 



1.0 



FIG. 9. The rotationally induced changes in the neutron number density, 770 and r/2, for our model star and the case 
when the superfluid neutrons corotate with the protons ($!„ = f2 p ).This particular data corresponds to the protons rotating at 
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FIG. 10. The rotationally induced change in the proton number density, $o and $2, for our model star and the case when 
the superfluid neutrons corotate with the protons (fi n = fl p ). This particular data corresponds to the protons rotating at 



Up = 1 kHz. 



Once we have solved the I = 2 equations and determined the value of A we can calculate the quadrupole moment 
from (ml). In doing this it is useful to construct two dimensionless quantities 



\GM J M 



and 



_ c J 

X ~ GAP 



(73) 



Once we do this, we can compare our results to those of Laarrakers and Poisson |3q | who study the quadrupole moment 
of rapidly rotating single fluid stars for realistic equations of state. In figure |11| we compare our results to theirs by 
showing the dependence of the ratio q/x 2 on the mass of the star. Clearly, the present results for the quadrupole 
moment are reasonable. It should also be noted that one would have q/x 2 = ~ 1 f° r a Kerr black hole. Thus figure fir] 
emphasizes the well-know fact that rotating neutron stars have a multipole structure that is radically different from 
that of the Kerr spacetime. 
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FIG. 11. Our results for the quadrupole moment (in terms of the ratio q/x 2 , see the main text for definitions) as a function 
of the stellar mass are compared to numerical results obtained for rapidly rotating relativistic stars using realistic supranuclear 
equations of state. 
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V. SEQUENCES OF STARS 

Given the machinery we have developed in the previous sections we are set to ask physically relevant questions. 
Ideally, we would like to be able to construct a sequence of rotating models that can meaningfully be said to represent 
the same non-rotating star. Given that we can allow the neutrons and the protons to rotate differently, and also 
adjust the relative central densities, the construction of such a sequence of models is non-trivial. However, a natural 
starting point may be to first determine the maximum allowed rate of rotation. Then we will discuss conservation 
of the total integrated neutron and proton numbers along a sequence of rotating stars. To impose this requirement 
corresponds to determining the total baryon numbers for the neutrons and the protons and then finding models such 
that these quantities are held constant as the spin-rate increases. Finally, we will turn to issues regarding chemical 
equilibrium. 

A. The Kepler limit 

The maximum spin-rate of a stationary configuration corresponds to the so-called Kepler frequency. It is determined 
as the rotation rate where mass-shedding sets in at the equator. The corresponding frequency can be calculated by 
requiring that the rotation rate of the fluid be equal to the frequency of a particle orbiting the star on an equatorial 
geodesic. The determination of the Kepler frequency for our slowly rotating stars is relatively straightforward. After 
all, the main part of the calculation involves the rotation frequency of an orbiting particle, and as far as the exterior 
vacuum is concerned, the fact that we are dealing with the rotation of two coupled fluids is irrelevant. The only 
subtlety of the two-fluid problem involves the fact that we can allow the protons and neutrons to rotate at different 
rates. Clearly, the Kepler limit will be determined by whichever fluid is rotating the fastest at the surface. In other 
words, we simply take Qk — max (fi n , ftp) in the following. It should, of course, be noted that the situation is simpler 
for more realistic neutron star models. Then the Kepler frequency will be determined from the rotation rate of the 
crust. 

For a rapidly rotating star the Kepler frequency VLk can be determined from (cf. Friedman, Ipser and Parker [p9[) 

Nv 
Vl K = ^=+u (74) 



K 



where the metric variables are defined by (pf) and 



K 3 ' 2 J 2KN> (K 3 ! 2 



NK' V NK' V NK 



*W 



(75) 



is the orbital velocity according to a zero-angular momentum observer at the equator (where all quantities should 
be evaluated). Here primes denote derivatives with respect to f. In order to be consistent, we expand this equation 
keeping only terms up to (and including) order Q 2 . Then we get 



, , h' fk' f 3 (u/) 2 

h-k + — ' 



v' 2 4v'e" 



0(n 3 ) . (76) 



The next step corresponds to using the vacuum expressions for the various slow-rotation quantities, and also doing 
the coordinate transformation f — v r. Luckily, the latter affects only the first term in (ffq). After some algebra we 
arrive at our final expression 

[Mm [MJSM (R + 3M)(3R~2M) p 3 2Jq + 6 *\ o2 {77 , 

V R 3 R s V R \ 2M 4R 4 M 2 4 R 2 J p V ' 

where we have made the scaling with Q p explicit by introducing J = JQ p etc. We have also defined 

3(i? 3 -2A/ 3 ), / 2M\ 3R 4 -3R 3 M -2R 2 M 2 - 8RM 3 + 6M 4 , s 

■kg U--5- + KWTTzrB o»^ ' ( 78 ) 



AM 3 °V R J 2RM 2 (R-2M) 

Here it should be noted that, at the ft 2 level the Kepler limit depends explicitly on the rotationally induced change 
in mass, the centrifugal flattening as well as the change in quadrupole moment (through A). 
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In order to determine the Kepler limit for our rotating superfluid models we simply have to solve the above quadratic 
for fl p . This then tells us the maximum allowed spin-rate for the protons at any given central energy density and 
relative (constant) rotation rate J7 n /fi p . Typical results of this exercise are shown in figure 12. In this figure we show 
how the Kepler frequency is affected by changes in the relative rotation rate f2„/J7 p . The results can be understood 
in the following way: For fl n /fl p = 1 the two fluids rotate together so we have a fairly standard result. For the model 
star from figure |^ we find that the Kepler frequency is Qk — ^l p — Sl n ~ 1-02 x 10 4 rad/s or vk ~ 1630 Hz. When we 
compare this to the relation (£2|) we see that it is not an unreasonable result, although it is roughly 10% higher than 
the empirical estimate for the relevant mass and radius. Figure |l2j also provides information regarding the effects of 
differential rotation. First of all it is clear that when £l n > tt p the Kepler rate is determined by the neutrons. Given 
that the neutrons make up more than 90% of the material in the star it is not too surprising that the result for £Ik 
changes very little as we increase O n beyond f2 p . Basically, the star rotates like a single fluid star with a dynamically 
small proton component. The opposite case, when f2„ < Q p , is more interesting. Then we see that the maximum 
allowed spin rate for the protons increases as Q n decreases. We can deduce that this is reasonable from (J77|): As we 
decrease the rotation rate of the neutrons the various VL p terms all decrease (again, because the neutrons make up 
most of the star and so have a dominant effect). Thus it follows naturally that the orbital frequency of a particle 
at the equator will approach the non-rotating star case. In other words, it increases which leads to the maximum 
allowed rotation of the protons increasing as well. 
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FIG. 12. The Kepler limit fix is shown as a function of the relative rotation rate fi n /fi p for our model star. The filled 
squares show the maximum allowed rotation rate for the protons (fi p ), while the open squares show the corresponding neutron 
spin rate (f2 n ). The Kepler frequency simply corresponds to the largest of the two. 



B. The baryon numbers and angular momenta 



We want to determine the conserved neutron number N n , 
total baryon mass Mb as 



the conserved proton number N p , and then define the 



M B = m n N n + m v N v 



(79) 



where m n is the neutron mass and m p is the proton mass. One of the assumptions behind the general relativistic 
superfluid formalism implemented here is that the neutrons and the protons are independently conserved. Thus the 
formulas for the conserved particle numbers can be obtained by integrating the conservation rules ( |ll|) over a bounded 
region of spacetime, and then using the generalized Stokes Theorem to reduce the spacetime volume integral to three- 
dimensional integration over the boundary. The bounded spacetime region is the four-dimensional volume contained 
between the two spacelike slices given by t = t\ and t = £2 where t\ and £2 are constants. The timelike boundary of 
this volume is the union of each sphere at infinity that exists on each t = constant spacelike slice. Because the matter 
is assumed to have compact support, there will be no contributions to the baryon mass formula from the timclike 
boundary. All that remains are two integrals, one on the t = ti and another on the t = £2 hypersurface, that add to 
zero. The t — ti integral, say, then yields the conserved particle number. 
The conserved neutron particle number is thus found to be 
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N n = 



n'*ri li yfyd 3 x 



(80) 



where 7 = det7y and 7^ is the intrinsic metric of the t = constant hypersurface, i.e. the metric is assumed to have a 
3 + 1 decomposition of the form 



ds 2 = - (N 2 - jijN 1 ^) dt 2 + 2 llj N j Ax l dt + 7ydar i cto j , 



(81) 



and the (future-directed) normal to the t = constant hypersurfaces, to be denoted ry M , is given by r/ M = (—AT, 0,0, 0). 
The formula for N p is identical except that n^ is replaced with j/\ One important point to keep in mind is that the 
surface of the star (which will be non-spherical in general) is given by r — constant, cf. Figure ffl. In terms of the 
metric and particle number density currents developed for the slow rotation formalism it is found that 



N n = 4n drr 2 e x/2 n ( 1 + tj + v Q + 



~2 



£0 



2r 2 f: 



(82) 



N p = 4?r 



I drr 2 e x/2 p ( 1 + $ + v + 



2 



2r 2 ~ 2 

e ° + 3^ 



(83) 



As an illustration we note that in the case of the model in Figure 0, and f2 n = J7 p , the two baryon masses are 

M Bn /M Q » 1.46 + 0.13 (j^Y 
M Bp /M q k 0.11 + 0.027 (j^Y ■ 

In the standard single fluid case one would typically identify a distinct sequence of rotating models as representing 
"the same" star by requiring that the total baryon rest mass be conserved. In our present problem, the fact that we are 
dealing with two distinct fluids provides additional complications. Provided that the star spins up (or down) rapidly 
enough compared to the nuclear reactions that convert neutrons to protons (and vice versa) it would be reasonable 
to conserve the total integrated neutron and proton neutron numbers individually. However, the spin-evolution of 
a neutron star typically proceeds slowly. Furthermore, as discussed by Reisenegger jl(j, one would expect a certain 
amount of conversion between the two fluids as the star spins down simply because of the change in shape of the 
constant density surfaces. In other words, in the two fluid case it is not obvious how to construct an "astrophysical 
sequence" of rotating stars based on the two baryon masses. 

Anyway, let us suppose that we want to demand that the individual integrated baryon numbers are conserved as 
we spin the star up. From the above results we see that this corresponds to 



5N n = / drr 2 e x/2 n Q ( m + v Q + 



A' 2 

-2 + ~r 



2r 2 - 



(84) 



and 



SNr, 



f drr 2 e x / 2 Po U 



■vo 



A' 2 



t 2r 2 - 2 
e ° + 3P L P 



(85) 



Let us consider the (by now familiar) stellar model from figure |2| According to the result above we then have 

M B = M Bn + Mb p = (1.46 + 0.11)M Q = 1.57M© , 

in the rotating case. First of all we find that the two baryon masses are not individually conserved for the star that 
has the same total M B at the Kepler limit. At least not if we take ?7o(0) = and Vt n /Vl p — 1. For this case we find 



M Bn = 1.48M, 



) 



and M Bp = 0.09M G 



However, we have the freedom to choose our parameters differently and we find that we can conserve both baryon 
masses by adjusting either i]q(0) or £l n /fl p (or, indeed, both of them). For example, we find that both M Bn and M Bp 
are the same as in the non-rotating star if we take 
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770(0) = and n n /n p = 0.75 



and Cl n /Q p = 1 



This is an interesting illustration of the subtleties involved in constructing rotating superfluid neutron star configura- 
tions. 

Furthermore, the fact that the slow-rotation calculation requires, once a non-rotating star is given, the specification 
of both Q n /Jlp and 770(0) leads us to ask an intriguing question: We have two seemingly free parameters and two 
equations, (p4) and (pq). Is it possible to specify these constants in such a way that the two baryon numbers are 
conserved even though we keep the central density of the star fixed? In the standard one-fluid case this would certainly 
not be possibly. As the spin of the star increases the centrifugal flattening will necessarily lead to decrease in the 
central density (for a star with a fixed number of baryons). Is it possible that two- fluid models offer a somewhat 
counter-intuitive alternative? That is, is it possible to carefully redistribute the neutrons and the protons (in such a 
way that the central energy density remains unchanged), and spin the star up without causing the familiar increase 
in baryon mass? For our particular model equation of state we find a negative answer to this question, cf. Figure |l3| . 
However, it is an intriguing possibility and we cannot at this point argue why it could not happen for some other, 
perhaps more realistic, equation of state. 
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FIG. 13. Curves corresponding to conservation of the total number of neutrons (SN n = 0) and the total number of protons 
(SN P = 0) are shown as functions of the two free parameters ?7o(0) and Q„/Q p . (The data corresponds to the case v v =1 KHz.) 
A point of intersection would correspond to a rotating model which conserves both the individual baryon masses from the 
non-rotating star. As is clear from the data, such models cannot be constructed for our model equation of state (which agrees 
with the standard result from the one-fluid case). This does not, however, rule out this possibility for other equations of state. 



Using the above set-up for the bounded spacetime region, we can also make meaningful definitions for the individual 
neutron and proton total angular momenta, to be denoted J n and J p respectively. Langlois et al (q| demonstrate 
that such an unambiguous separation can be made if the hypersurface over which the integrals are being performed 
is invariant under the axisymmetry action. In the present context this means that the scalar product of the spacelike 
Killing vector ^ with the normal to our hypersurface 77^ must be zero. Indeed, it is trivial to verify that we have 
(\> v r\ v = 0. Now, according to Langlois et al, the two angular momenta are given by 

(86) 



(87) 



J t— const J t— const 

In terms of the slow-rotation approximation, we find for the neutron total angular momentum 



J n = -5E f drrV*-)/ 2 
<J Jo 



and the very similar form 



87T 



drr 4 e (A-„)/2 



fj, n L/ n -j- ^.qTIoPo \* ifi ^Zp) 



XoPoL p + A n p (fl p - fi n ) 
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for the proton total angular momentum. It is easy to verify that 

J = J n + J p , (89) 

where J is the total angular momentum found earlier by matching to the exterior vacuum solution of Hartle. 

C. Chemical Equilibrium 

So far we have not imposed chemical equilibrium in our rotating stars (although it was assumed for the non-rotating 
background model) . This is quite reasonable since one can easily think of situations where the two fluids are temporarly 
rotation "out of sync" following for example a glitch. As we argued in the Introduction it will then likely take the 
various relevant mechanisms many hundreds of dynamical timcscalcs to again lock the two components together. And 
it may take much longer for chemical equilibrium to be restored given that this happens on the characteristic timescale 
of the various nuclear reactions. Anyway, before we conclude this paper it is obviously meaningful to address the issue 
of chemical equilibrium. Intuitively, one might expect that one can only have equilibrium when the neutrons and the 
protons corotate. The argument for this is simple: If the two species rotate at different rates we must work in one of 
the two frames (the nuclear reactions that re-establish equilibrium depend on the local physics), and the result will 
depend on which frame we choose. In the frame rotating with the neutrons the proton chemical potential will have 
an additional kinematic piece. This is also true for the neutron chemical potential in the frame that corotates with 
the protons. In other words, it would seem possible to have equilibrium only if the two fluids corotate. 

This is, of course, just a hand-waving argument and we need to support it mathematically. Let us first return to a 
situation where fi„ and fi p are not assumed to be constant. Langlois et al jq] have shown that chemical equilibrium 
is imposed on the system via the condition 

f ijh, - Xv) = ■ (90) 

Using the results from Section II we immediately find that this corresponds to 

t v n» + Slpffiv = t v X v + Sl v <\> v Xv ■ (91) 

Let us now assume, again, that fi p is constant. As we have already argued in Section IV, this is likely to be the case 
in most astrophysical neutron stars. The Euler equation for the protons (n2h then implies that the right-hand-side of 
the above condition is a constant. Thus we have 

t v Xu + &p4> v Xu = constant . (92) 

Obviously, this also implies that 

t v n v + Qp4> v ii v = constant . (93) 

We solve for t v \i v and then put the result into the Euler equation for the neutrons. This way we obtain 

(fin - «p) Wa**) = . (94) 

From this we can conclude the following: If one demands chemical equilibrium, and also assumes that fi p is constant 
(that the protons rotate rigidly), then one of the following two conditions must be true, 

fi„ - fi p = , (95) 



4> v Hv = constant . (96) 

The first of these is the result that we anticipated, that chemical equilibrium would lead to the two fluids co-rotating. 
The second condition is a bit more puzzling, and we need to ask whether one can have a physical configuration that 
satisfies it. 

Taking the slow-rotation form for \i v the second condition becomes 

- r 2 sm 2 6e~ v/2 [B n o {u - fi„) + A p {u - fip)] = constant . (97) 
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From this it is clear that unless the right-hand-side vanishes then u> and f2„ are both badly behaved at the origin (as 
well as at the poles). Thus we can only have a well-behaved solution if 

Q n = [(B n + A oPo )/(B n )}uj - [A oPo /B n ]n p . (98) 

In the special case with no entrainment we have A = 0, and our condition implies 

Q n = u> . (99) 

Hence, we see that the second condition for chemical equilibrium condition links the rotation rate of the superfiuid 
neutrons to the frame-dragging. This means that, in order to satisfy this condition, the neutrons must rotate differ- 
entially. At first this possible solution may seem peculiar, but it may in fact make sense if we consider the physical 
meaning of the frame-dragging. The condition (|99| ) simply says that the neutrons are not rotating with respect to a 
local zero-angular momentum observer. They are simply dragged along by the rotation of the protons. 

From this discussion we conclude that if the two fluids both rotate uniformly and tt n ^ tt p then one cannot impose 
chemical equilibrium. When chemical equilibrium is imposed, the two fluids must either corotate, or fi„ must have 
differential rotation (and in fact be equal to the frame-dragging m when there is no entrainment). 

We note that one could, in principle, explore the latter possibility further. We can combine the condition for 
equilibrium with the frame-dragging equation (|4j) in an interesting way. First we re-express (|44|) in terms of L p . 
Then we get 

^ (rV^/ 2 ^)' - I6^ A -^/ 2 (tf - A ) L p = I67re^/ 2 ti n (fi p - Q n ) . (100) 

By making use of (J9q) we readily rewrite this as 



1_ Og-CA+.O/a^y _ le^CA-O/a 



*o - A + ^ (B n + A oPo ) 



L p = . (101) 



This alternative frame-dragging equation now takes the same form as the single-fluid equation derived by Hartle |3C|] 
(the difference being in the factor multiplying L p ). Thus we know that if we were to solve it numerically we would 
find that L p decreases monotonically from the centre of the star to the surface. In other words, the results would be 
similar to those shown in figure pi In view of this we do not show such results here. 

VI. CONCLUSIONS 

In this paper we have developed a framework for constructing and analyzing slowly rotating relativistic superfiuid 
neutron stars. The model is based on the standard two-fluid representation for superfluids wherein all the charged 
components: protons, electrons and lattice nuclei, are considered as a single fluid coexisting with the superfiuid 
neutrons. We have applied the formalism to a simple equation of state for which the two fluids are described by 
polytropes. We have then studied the effects of rotation on the resultant stars, with particular focus on the effects 
due to the fact that the two fluids need not rotate together. 

The present results provide a framework that opens the door to fully relativistic studies of many important astro- 
physical problems. We are currently extending this work in two directions. First, we are investigating the effects of 
entrainment on rotating stellar models. Secondly, we are extending previous studies of quasinormal mode oscillations 
in relativistic superfiuid stars p8[ to incorporate the effects of slow rotation. Of particular interest would be a study 
of the low-frequency inertial modes, following in the footsteps of Lockitch et al. [g6|. The calculation of such modes 
has been brought into sharp focus since the discovery that the r-modes are generically unstable due to the emission of 
gravitational waves. So far the only study of the r-modes in the superfiuid context is the Newtonian work of Lindblom 
and Mendell p8fl . Our plan is to study the same problem within the framework of relativity. This is an issue that 
demands immediate attention given the possibility that the r-modes may lead to detectable gravitational waves and 
the fact that a new generation of large scale interferometric detectors will come on-line soon. 
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VII. APPENDIX I 

Here we give the four combinations used in the right-hand-sides of the field equations that determine the four metric 
cocffcients 
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Some relationships useful for obtaining the stress-energy combinations just given are the following: 
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with similar results for B and C, 
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(104) 
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(105) 



and 
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where ^ — A = /j, n + XoPo- The analogous expressions for the Einstein and Ricci tensor components can be found 
in Hartle iotl. 
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